clear all
%close all
addpath('sub_Table')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Setting: Parameters
% Generating artificial data from non-homothetic preference
paras.flag_sigma = 0; % 0 for constant sigma, 1 for variable sigma
N = 5000;
I = 3;% Need to adjust it when the number of goods changed.
T = 100;
pols = [1,2,4,6,8,12];

name = 'const_support';
% Nh-CES parameters
if paras.flag_sigma == 0
    paras.sig = 5;
else
    paras.sig = 10;
end
paras.eta = -2;
paras.eps_vec = [0.8;1;1.65]*2;%
paras.T = T;paras.T = T;paras.N = N;paras.I = I;%paras.Ns = Ns;

%% Setting: Price and Wage
% Price Schedule
pvec(1,1,:) = reshape(exp(linspace(0,log(2),T)),[1,1,T]);% (I,1,T)
pvec(2,1,:) = reshape(exp(linspace(0,log(3),T)),[1,1,T]);% (I,1,T)
pvec(3,1,:) = reshape(exp(linspace(0,log(4),T)),[1,1,T]);% (I,1,T)


% Generate artificial util:  
trendU1 = (reshape(exp(linspace(0.3,2,T/2)),[1,1,T/2]));
trendU2 = (reshape(exp(linspace(2,1,T/4)),[1,1,T/4]));
%trendU3 = (reshape(exp(linspace(1,0.3,T/4)),[1,1,T/4]));
trendU3 = (reshape(exp(linspace(1,0.25,T/4)),[1,1,T/4]));


trendU = cat(3,trendU1,trendU2,trendU3);

trendL1 = (reshape(exp(linspace(0,-2,T/2)),[1,1,T/2]));
trendL2 = (reshape(exp(linspace(-2,-0.5,T/4)),[1,1,T/4]));
trendL3 = (reshape(exp(linspace(-0.5,0.2,T/4)),[1,1,T/4]));

trendL = cat(3,trendL1,trendL2,trendL3);

trendVec = zeros(1,N,T);
for t = 1:T
    trendVec(1,:,t) = linspace(trendL(:,:,t),trendU(:,:,t),N);
end
v_vec_tmp = linspace(0,1,N).*(trendU- trendL) + trendL;

v_vec = sort(v_vec_tmp,2);

paras.omega = [1;1;1];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Generating Artificial data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[I_vec, u_vec, B_vec]= NhCESmaxV(v_vec,pvec,paras) ;

%%BBK Algorithm %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[U_vec_FP] = CalMoneyMetric(I_vec, B_vec, pvec,1);
[U_vec] = CalMoneyMetric(I_vec, B_vec, pvec,0);
%% Run JL with loop
format shortG

for i = 1:size(pols,2)
pol= pols(i);    
paras.Kn = pol;
warning("off")
 lastwarn('');
[LJMM_F,LJMM_SC,JLerrorType] = CalLJ_poly(I_vec, B_vec, pvec,paras);
run("sub_err.m")
ErrorJLVec(i,:) = ErrorJL;
JLErrorTypeVec(i,:) = JLerrorType;
end

%% produce tables
caption_text = 'JL Algorithm';
latex_table = create_latex_table_JL(ErrorJLVec,pols, JLErrorTypeVec, caption_text)
%save_latex_table_to_file(latex_table, ['../fig/' name '_JL.tex']);
caption_text = 'BBK Algorithm';
latex_table = create_latex_table_1x4(ErrorBBK);
%save_latex_table_to_file(latex_table, ['../fig/' name '_BBK.tex']);


function [I_vec, u_vec, B_vec]= NhCESmaxV(v_vec,pvec,paras) 
%%Nonhomothetic Preference: Given Wfull and P, solve for u 
T = paras.T;
N = paras.N;
I = paras.I;
%Ns = paras.Ns;

I_vec = NHU(v_vec,pvec,paras) ;
%%Obtain EV
[u_vec,B_vec]= CalEV(v_vec,pvec,paras);

function I_vec = NHU(v_vec,pvec,paras) 
if paras.flag_sigma ==1
    sig0 = paras.sig;
    xi = 0.01;
    xi2= 100;
    paras.sig = ((1.5)^((xi2-1))+((5^((xi-1)/xi)+(10-2*log(v_vec)).^((xi-1)/xi)).^(xi/(xi-1))).^((xi2-1))).^(1/(xi2-1));
end
I_vec = sum(paras.omega.*(((v_vec.^paras.eps_vec).*pvec).^(1-paras.sig))).^(1./(1-paras.sig));
end

function [u_vec,B_vec]= CalEV(v_vec,pvec,paras)
paras.sig =paras.sig;
sig0 = paras.sig;
if paras.flag_sigma ==1
    xi = 0.01;
    xi2= 100;
    paras.sig = ((1.5)^((xi2-1))+((5^((xi-1)/xi)+(10-2*log(v_vec)).^((xi-1)/xi)).^(xi/(xi-1))).^((xi2-1))).^(1/(xi2-1));
end
pvec_b = pvec(:,:,1);
u_vec = sum(paras.omega.*(((v_vec.^paras.eps_vec).*pvec_b).^(1-paras.sig))).^(1./(1-paras.sig));
tmp = paras.omega.*(((v_vec.^paras.eps_vec).*pvec).^(1-paras.sig));
B_vec = tmp./sum(tmp);
end

end